function ImageBinary=PlotSynergy(L4,L5,epoch,FN_i,CameraAngleCapability,DU,VU,TU)

options=odeset('RelTol',1e-13,'AbsTol',1e-13);


%% Compute L4

S_i = L4.S_i;

temp_time_month = month(L4.epoch);

ind = find(temp_time_month==epoch);

% plot3(S(1,:),S(2,:),S(3,:))
[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
parfor ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility_i(ii,jj) = 1;
            end
        end
    end
end


ImageBinary_L4 = mat2gray(visibility_i);

Window = figure('OuterPosition',[0,0,500,950]);
subplot(3,1,1); hold on;
mesh(ImageBinary_L4*100,'EdgeColor','none','FaceColor','flat');
xlabel('Longitude')
ylabel('Latitude')
title('L4')
set(gca,'fontsize',16)
xticks(linspace(0,row,5));
xticklabels({'180W','90W','0','90E','180E'});
yticks(linspace(0,col,5))
yticklabels({'90S','45S','0','45N','90N'});
xlim([1,row]);
ylim([1,col])

%% Compute L5

S_i = L5.S_i;

temp_time_month = month(L5.epoch);

ind = find(temp_time_month==epoch);

% plot3(S(1,:),S(2,:),S(3,:))
[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
parfor ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility_i(ii,jj) = 1;
            end
        end
    end
end


ImageBinary_L5 = mat2gray(visibility_i);

subplot(3,1,2);hold on;
mesh(ImageBinary_L5*100,'EdgeColor','none','FaceColor','flat');
xlabel('Longitude')
ylabel('Latitude')
title('L5')
set(gca,'fontsize',16)
xticks(linspace(0,row,5));
xticklabels({'180W','90W','0','90E','180E'});
yticks(linspace(0,col,5))
yticklabels({'90S','45S','0','45N','90N'});
xlim([1,row]);
ylim([1,col])
%% L4+L5

ImageBinary_L45 = ImageBinary_L5+ImageBinary_L4;

subplot(3,1,3);hold on;
mesh(ImageBinary_L45*100,'EdgeColor','none','FaceColor','flat');
xlabel('Longitude')
ylabel('Latitude')
title('L4+L5')
set(gca,'fontsize',16)
xticks(linspace(0,row,5));
xticklabels({'180W','90W','0','90E','180E'});
yticks(linspace(0,col,5))
yticklabels({'90S','45S','0','45N','90N'});
xlim([1,row]);
ylim([1,col])




end